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Abstract 

The benchmark active controls technology and 
wind tunnel test program at NASA Langley Research 
Center was started with the objective to investigate the 
nonlinear, unsteady aerodynamics and active flutter 
suppression of wings in transonic flow. The paper will 
present the flutter suppression control law design process, 
numerical nonlinear simulation and wind tunnel test 
results for the NACA 0012 benchmark active control 
wing model. The flutter suppression control law design 
processes using classical, and minimax techniques are 
described. A unified general formulation and solution for 
the minimax approach, based on the steady state 
differential game theory is presented. Design 
considerations for improving the control law robustness 
and digital implementation are outlined. It was shown that 
simple control laws when properly designed based on 
physical principles, can suppress flutter with limited 
control power even in the presence of transonic shocks 
and flow separation. In wind tunnel tests in air and heavy 
gas medium, the closed-loop flutter dynamic pressure was 
increased to the tunnel upper limit of 200 psf. The control 
law robustness and performance predictions were verified 
in highly nonlinear flow conditions, gain and phase 
perturbations, and spoiler deployment. A non-design 
plunge instability condition was also successfully 
suppressed. 

Introduction 

The benchmark active controls technology 
(BACT) and wind tunnel test program at NASA Langley 
Research Center was started with the objective to 
investigate the nonlinear, unsteady aerodynamics and 
active flutter suppression of wings in transonic flow. 
Under the initial wind tunnel test program, a NACA 0012 
airfoil rectangular wing, equipped with pressure 
transducers, active trailing edge control surface, and two 
spoilers were constructed for active flutter suppression 
tests. The model was mounted on a pitch and plunge 
apparatus in the NASA Transonic Dynamics Tunnel in 
order to test flutter suppression control laws and measure 
unsteady pressure distributions in nonlinear flows with 
oscillating shocks and boundary layer separation. It was 
necessary to develop a flutter suppression system that 
would be stable under these flow uncertainties. 


This paper describes flutter suppression control 
law design processes using classical and unified linear- 
quadratic Gaussian minimax techniques. A unified 
general formulation for the linear quadratic Gaussian and 
minimax methods based on the steady state differential 
game theory is presented. Lessons learned in evaluating 
and improving the singular value based multi-input multi- 
output system robustness are described. Design 
considerations for digital implementation are outlined. 
Numerical simulation of the control law performance, and 
wind-tunnel test results for flutter suppression, are also 
presented. 



Fig. 1 BACT model test setup in wind tunnel 

Wind Tunnel Model Description 

A perspective view of the BACT model test set 
up on the Pitch and Plunge Apparatus (PAPA) in the wind 
tunnel is shown in Fig. 1 . Fig. 2 shows the control surface 
and sensor locations. The rigid wing section has pitch and 
plunge degrees of freedom. The accelerometer sensors are 
located near the section leading edge ( zle ) and trailing 
edge (zte) at the section inboard. An identical pair of 
sensors is located at the section outboard as a spare. The 
partial span spoilers are located on the upper and lower 
surfaces, just ahead of the trailing edge control surface. 
Each of the control surfaces stretched over 30% of the 
span and 25% of the chord. The bending and torsion 
frequencies of the PAPA mounted NACA 0012 wing 
model were 3.3 Hz and 5.2 Hz respectively. 
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diverges at the rate of 6 lbs/sec. The moment diverges at a 
rate of 1 lb/sec. 


(Pitch and Plunge Apparatus) 



Fig. 2 NACA 0012 BACT wing on PAPA. 

Preliminary Analysis 

The preliminary analysis, control surface sizing, 
and flutter suppression control law design were based on 
the analytical state-space equations of motion of the 
BACT wing model. 1 4 These equations were developed 
analytically, using structural dynamic analysis and 
unsteady doublet lattice aerodynamics with rational 
polynomial approximations 5 . These linear state space 
equations consisted of 14 States (plunge, pitch, plunge 
rate, pitch rate, 3 aerodynamic states for plunge, 3 
aerodynamic states for pitch, 2 trailing edge flap actuator 
states, 2 Dryden gust states), 2 inputs (actuator command 
and gust input noise) and 7 outputs (zte and zle 
acceleration, flap command, flap deflection, rate, 
acceleration, and gust velocity). This 14 th order state 
space equation was used for classical control law design 
and for performance simulation and verification purposes. 
For the optimal control law design purposes and for 
presentation of the design data in a concise form, the 14 th 
order state-space equations were reduced to 4 th state-space 
equations, using residualization and Schur's balanced 
reduction method 6 7 . First, it was reduced to an 8th order 
system using residualization technique, in which only the 
static part of all modes above 15 Hz were retained. The 
resulting 8th order system was then balanced and the four 
states of the system with largest balanced singular values 
were retained. A sample of the 4 th order model design data 
is presented in the Appendix. 


Open loop TE Accl. 225 psf Open Loop LE Accl. 225 psf 



Time, sec Time, sec 

Fig. 3 Open loop transient responses in air at 225 psf. 
Frequency Responses 

The open loop frequency responses were studied 
using this 14 th order plant model, to select a possible 
candidate for feedback signal in the flutter suppression 
control law design. The Bode diagram of the trailing edge 
and leading edge accelerometers (zte) and (zle) and their 
difference (zte- zle) due to the trailing edge control 
surface excitation (6te) in air at 225 psf at Mach 0.5, are 
shown in Fig. 4. The magnitude plots indicate 
predominant plunge response at 3.3 Hz excitation 
frequency. At 4.2 Hz excitation, the motion is a 
combination of pitch and plunge with pitch motion 
leading the plunge. The (zte -zle) represents a signal 
proportional to the pitch acceleration and can be 
integrated to provide a pitch-rate signal. Feedback of this 
signal with proper gain can provide maximum pitch 
damping at the flutter frequency. 


zte/dte, zle/dte, (zte-zle)/dte 225 psf 



Open-loop Responses 

The analytical open-loop flutter dynamic 
pressure in air was 128 pounds per square feet (psf) at a 
flutter frequency of 4.5 Hz. Fig. 3 shows the response of 
the wing trailing edge and leading edge accelerometers 
due to a 1 degree step input of the trailing edge control 
surface in air at 225 psf dynamic pressure. The primary 
plunge motion mixed with small pitch diverges rapidly. 
The unsteady lift forces oscillate about 8 lbs mean lift and 



Fig. 4 The Bode diagrams of zte and zle and (zte - zle) 
due to dte excitation in air at 225 psf, Mach 0.5. 
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Classical Control Law Design 

Based on this Bode plot, a classical flutter 
suppression scheme using pitch-rate proportional 
feedback from the zte and zle accelerometers was first 
devised by studying the Nyquist diagrams. The Nyquist 
diagram of the difference between trailing edge and 
leading edge accelerometers (zte - zle) due to the trailing 
edge control surface excitation (Ste) in air at 200 psf, is 
shown in Fig. 5(a). The arrow indicates increasing 
frequency of excitation from 2 FIz to 6 Hz, with each * 
representing frequency increment of 1 radian/second. 
Since the open-loop plant had a pair of complex unstable 
poles, and the Nyquist contour did not encircle the -1 
point, the unit feedback closed-loop system would be 
unstable. However, if the (zte - zle) signal was integrated 
to provide a 90 degree phase lag and then used for 
feedback with sufficient gain, the Nyquist contour would 
rotate 90 degrees clockwise and then expand to encircle 
the -1 point to achieve stability. A washout filter of type 
s/(s + a) was also required, to remove any static bias that 
would otherwise be amplified by the integration. The 
series connection of integrator 1/s and washout filter was 
equivalent to a first order lag filter a/(s + a), where s is 
the Laplace operator. 


zte, g 



(zte-zle) 


zle, g 



Fig. 5(a) Nyquist diagram of (zte-zle) due to Ste 
excitation in air at 200 psf, Mach 0.5 


Gain Selection 

Two types of lag filters, namely 5/(s + 5) and 
10/(s + 10) were examined. The latter was selected to 
achieve a higher phase margin at the plant input above the 
flutter frequency. Higher phase margin was desirable for 
two reasons 7 . First, the 25 Hz antialiasing filter and the 
1/200 seconds computational delay contribute about 20 
degrees of phase lag at the flutter frequency. Secondly, 
with increasing dynamic pressure, the actuators may have 
additional unknown phase lag, as the control surface 


moves against higher aerodynamic loads. The Nyquist 
diagram of the difference between trailing edge and 
leading edge accelerometers (zte - zle) with 10/(s + 10) 
lag filter and a gain KR = 500 due to the trailing edge 
control surface excitation (Ste) in air at 200 psf, Mach 0.5, 
is shown in Fig. 5(b). The unit circle is also shown. 
Because the Nyquist contour encircled the -1 point, the 
unit feedback closed loop system would be stable. As 
desired, the phase margin at the plant input above the 
flutter frequency was about 60 degrees, but the phase 
margin below the flutter frequency was only 20 degrees. 
Preliminary analysis indicated that this basic simple 
control law 1 


Ste = 500 — (zte - zle) 

s + 10 

can suppress the flutter instability in the dynamic pressure 
range from 0 to over 225 psf, both in air and in heavy gas 
medium. However, the closed loop transient responses 
and stability margins required substantial improvement. 


zte, g 



zle, g 



Fig. 5(b) Nyquist diagram of (zte -zle) with 10/(s+10) 
lag filter and a gain KR = 500, due to Ste excitation, in air 
at 200 psf, at Mach 0.5 . 

Root Locus 

Analysis of the root locus with pitch 
acceleration (zte - zle ) feedback through a 10/(s + 10) lag 
filter with increasing gain KR = 0, 500 , ..., 2500 is shown 
in Fig. 6(a). The stabilization was achieved by increasing 
the pitch model damping and lowering the plunge mode 
frequency. An additional feedback of the pitch rate 
proportional signal through a 5/(s + 5) lag filter with KR 
= 500 and increasing gain KP =0,500, .... 2500 was used 
to increase the damping and frequency separation further, 
as indicated by the root-locus diagram shown in Fig. 6(b). 
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This design strategy was equivalent to pitch-angle and 
pitch-rate proportional feedback that increased the pitch 
mode frequency and plunge mode damping. 


imaginary part, 
rad/sec 


imaginary part, 
rad/sec 



<5, =[500 1] 1 L 


rea l P ar t ra( t/sec rea l p ar ^ rac )/ sec 

Fig. 6(a) Root-locus with (zte - zle) feedback through a 
10/(s + 10) lag filter, with increasing gain KR (at left). 
6(b) Root-locus with additional pitch-rate feedback 
through a 5/(s + 5) lag filter with KR = 500 and 
increasing gain KP. The arrows indicate increasing gain. 

Pitch and Pitch-Rate Feedback Control Law 

From the root-locus study, the feedback gains 
were selected as KR = 500 and KP = 1. Thus, the second 
order state space equations of the initial pitch and pitch- 
rate feedback control law 2 is given by 


*ij r — 10 OlfXjj HO -lOlfe 

-*2 [ 5 -5\x 2 \ 0 0 }zf e 

x, 

X. 


The control law inputs are zte and zle in g unit 
and the output 8te is in degrees. The high feedback gain 
was required because the maximum (zte - zle) signal was 
only of the order 0.1 g/deg. The response exhibited 2% 
settling time in 1.5 seconds. However, the high gain 
resulted in a severe sensitivity with respect to plant 
perturbation and individual sensor uncertainty, as 
indicated by the corresponding singular value plots in Fig. 
7. Here G, K and A denote plant, controller and 
uncertainty block, respectively 8,9 . This figure indicates 
that the minimum singular value a(I+KG) is 0.3 at plant 
input and o(I+GK) is only 0.01 at plant output. This 
means that at 225 psf dynamic pressure, the closed-loop 
system has very little robustness to multiplicative 
perturbation 8 ' 9 at the plant output. 

These singular value a plots can be related to 
multivariable gain and phase margins using the universal 
gain and phase margin diagram 8 shown in Fig. 8. For 
example, minimum singular value a(I+KG) of 0.3 is 


equivalent to ±20 degrees phase and +3 dB gain margins 
at the plant inputs. The singular value l/o[K(I+GK) ’] is 
close to 0.005 g/degrees near 2 Hz. This means that the 
plant has very little tolerance to an additive perturbation A 
to the plant. The complex determinant locus of (I+KG) 
and its distance from the origin is a measure of its 
closeness to singularity. 



Fig. 7 Singular value plots for analysis of multivariable 
stability margins with a perturbation A, at the plant input 
or output, with classical control law 2, at 225 psf, in air. 


a Minimium singular value 



Fig. 8 Universal gain and phase margin analysis diagram. 

Final Pitch and Pitch-Rate Feedback Control Law 

This lack of robustness associated with this pitch 
and pitch rate feedback control law was alleviated by 
choosing a feedback of a proper linear combination of the 
two sensors with lower gains for KR, instead of using the 
difference (zte -zle). The linear combination of zte and 
zle, which is equivalent to feeding back both pitch 
acceleration (zte - zle) and plunge acceleration (zte + zle) 
in the ratio 0.7 (zte- zle) + 0.3(zte + zle), appeared to 
provide a superior control law. The final classical 
feedback control law, using this combination that is 
equivalent to (zte - 0.4 zle) feedback, along with reduced 
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gains of KR = 50 and KP = 1, was analyzed and 
implemented. The basic control is shown here in state- 
space form and is denoted by classical control law 3. 
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Response and Robustness Analysis 

The closed-loop transient responses due to 1 
degree step deflection of 8te, in air at 225 psf, at Mach 
0.5, is shown in Fig. 9. The trailing edge control surface 
shows only 0.25 degrees overshoot with a maximum rate 
of 12 degrees / sec. The lift and moment forces indicate 
about 20% load alleviation compared to the open-loop 
initial transient values shown in Fig. 3. Figure 10 shows 
the singular-value plots for analyzing the system stability 
margins 8,9 with law 3 at 225 psf dynamic pressure. Flere 
G, K and A denote plant, controller and uncertainty block 
transfer function, respectively. This figure indicates that 
the minimum singular value o(l+KG) is increased to 0.8 
at plant input and at plant output o(I+GK) is increased to 
0.3 from the corresponding values with law 2 presented in 
Fig. 7. The minimum singular value c(I+KG) of 0.8 is 
equivalent to +45 degrees phase and -5 dB to 12 dB gain 
margins at the plant inputs. These gain and phase margins 
are determined from Fig. 8 as previously described. The 
minimum singular value l/a[K(I+GK) ’] is also increased 
from 0.005 g/degree to 0.04 g/degree near flutter 
frequency, thus increasing the plant's tolerance to additive 
plant perturbation. The complex determinant loci of 
(I+KG) ideally should be outside the unit circle to achieve 
+6dB gain margins and +60 degrees phase margins. The 
computational delay and antialiasing filters added 20 
degrees phase lag. Hence, the system nearly attained these 
margins. The singular value plots indicate that the system 
is stable with adequate singular value based multivariable 
stability margins even at this high design dynamic 
pressure of 225 psf. This pressure is 97 psf above the 
open-loop flutter dynamic pressure q flutter 128 psf, 
representing a 75% increase. 


Unified Optimal Design 

Flutter suppression control law design using an 
unified (1) linear quadratic Gaussian (LQG) and (2) 
Minimax method 10 4l is presented next. The Minimax 
approach is analogous to the time domain H-infinity 
design 12 and is based on the steady state differential game 
formulation. The unified formulation of these optimal 
design techniques provide a basic understanding of the 
relation between them. The derivation from basic 
principles using variational principles are provided. The 


solution only requires an eigen-solver. The corresponding 
Matlab script is presented in the Appendix. 





Fig. 9 Closed-loop responses: control surface deflection 
and rate, lift and pitching moment, due to step input Ste 
with control law 3, at 225 psf, in air, at Mach 0.5 (open 
l0 °P V„er= 128 P Sf ^ 


-(l+A)-G-> -G-(l+A)-> 



Fig. 10 Singular value plots for analysis of multivariable 
stability margins to perturbation A, with classical control 
law 3, at 225 psf, in air. 


Unified Minimax Formulation 

Consider the state space Eqs.(l-3) representing 
the nth order plant, control input u(t), disturbance w(t), 
design output yd and sensor output y s , where all necessary 
rank, controllability and observability conditions are 
assumed to be satisfied. 


Plant state-space equations 
dx(t)/dt = F x(t) + G u(t) + G w w(t) 
and x(0) = x(j (1) 

Design output 

y'd(t) = H d x(t) + E du u(t) (2) 

Sensor output 

y s ( t) = H s x( t) + E sw w( t) (3) 
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State -Feedback Minimax Regulator Problem 

The Minimax problem is to determine the plant 
input u(t) which would minimize the quadratic 
performance index J , and find the worst plant 
disturbance w(t) and initial condition x„, which would 
maximize J defined in Eq. (4), 

°o 

J = if (x Q x x + x Q xu u + u Q u u)dt (4) 

subject to the constraint Eq.(l) with xq t xq = 1 and 
specified W defined by, 


Yt\ \ (F-G&Ql) (-GQ;'G T + r 2 G n K l Gl ) 1 J*(0 

' dX \ [-(Q x + Q xu Q: 1 Q t J -C f-gq;'q t J t J{A(0 

dt 

withx(0) = xo and Afco ) = 0 (12) 

State -Feedback Regulator 

Substituting A (t) = S(t)x(t), in Eqs. (10-12), leads 
to Eqs. (13-15). The general Riccati Equation (15) is then 
solved for the unknown n x n matrix S. 


u(t) = -Q u - 1 (G T S + Q xu T )x(t) 


w(t)= y~ 2 R w - ] G w T Sx(t) 


Usually, the constant weighting matrices Q x , Q u are 
unity, and Q xl{ = [0] in a H-infinity exposition. These are 
included herein to derive a unified general time-domain 
formulation. The cross weighting matrix Q xu originates if 
one uses y d from Eq. (2) in the performance index J to 
replace x. Then, Q x = H d T Q yd H d , Q xu = H d r Q yd E du , and 
Qu is replaced by [ Q u + E du r Q yd E du ]. The significance 
of the cross weighting matrix Q xu and how it can be 
selected for pole-placement of the state regulator will be 
shown later in the state-feedback regulator subsection. 

The minimax solution is given by the stationary 
condition of the augmented performance index J 

oo 

J = jf Q (x Q x x + 2x Q xu u + u Q u ii - y w R w w)dt 

oo f] Y 

+ A r J (Fx + Gu + G w w - —)dt + y 2 W (6) 

where, y is a scalar parameter. Using the calculus of 
variation with respect to x(t), u(t), w(t) and the vector 
Lagrange multiplier Aft), the conditions for dJ=0 are 
given by Eq.(l) and Eqs. (7) to (9). 


dh/dt 

- ~ Qx x ~ F^X - Q xu u 

(7) 

Qu lt 

- ~ G T X - Q x Jx 

(8) 


dS/dt + SF + F t S + Q x -(SG + Q xu )Qu- ] (SG + Q XU ) T 
+ S(y~ 2 G w R w ~ 1 G w r )S =0 (15) 

The positive definite symmetric solution for S is obtained 
from the (2n x n) eigenvectors of the n stable eigenvalues 
of the Flamiltonian matrix inside the square bracket [ ] in 
Eq.(12). For the steady state problem (i.e. dS/dt - 0 ), 
only the steady part of the Riccati Equation (15) is solved 
in order to obtain the symmetric positive-definite matrix 
S. If the eigenvectors are partitioned into two n x n 
matrices X and A which represent the stable subspace 
eigenvectors of x and A, then S = X~ J A. The constant 
optimal feedback gains, C 0 and C w , and the closed-loop 
system matrix are given by, 

C 0 =-Qu~ 1 (G T S + Q xu T ) (16) 

C w = y- 2 R w ~ 1 G w r S (17) 

dx/dt = [F + G W C W + GC 0 ]x. (18) 

Using Eqs. (1), (15) and (18), it can be shown 11 that 
optimal J and W defined in Eqs. (4) and (5) are given by, 

J =0.5 Trace [S] (19) 


y2R w w = G w T . A 


(9) W = 0.5 Trace [C W T R W C W X], 


(20) 


Solving for u(t) and w(t) from Eqs. (8-9) and substituting 
them in Eqs. (1) and (7), the necessary stationary 
conditions for J are obtained as. 


u =-Q u ~ 1 (G T X + Q xu r x) 

(TO) 

w = y~ 2 R w - J G w t X 

(11) 


where X is the solution of the Lyapunov Equation (21) 

[F+G W C W + GC 0 ] X + X[F + G W C W + GC 0 ] T + 
x(0)x(0) T = [0] (21) 

The worst x(0) that maximizes J is given by the 
eigenvector of the maximum eigenvalue of S. The 
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standard linear quadratic regulator (LQR) solution is 
obtained when 7 = 00 , (i.e. C w = 0). As y 2 is decreased, 
the worst response due to the disturbance w(t), measured 
by the maximum singular value of [x T Q / /2 u T Q u 1/2 [, is 
reduced. The minimum value of y 2 for which a stable 
solution of Eq. (15) exists provides the minimax state- 
feedback regulator that minimizes the maximum singular 
value of [x t Q x 112 u t Q u 1/2 ]. 

The State-Estimator Equation 

The derivation of coupled state-estimator 
equations using linear quadratic minimax approach is still 
a subject of research. Here the equivalent state-space 
solutions 12 of the H-infinity problem are presented. The 
state estimator gain B o D 0 is obtained by finding the 
symmetric positive definite solution for P from the state 
estimator Riccati Eq. (24) which is dual to the state 
regulator Riccati Eq. (15). 


G w H/ 

Qx G W R W G W 7 

Qu Rv 

s P 

Qxu Rwv 

(tf-t) (t-to) 

Table 1. Duality relations between linear quadratic state- 
regulator and state-estimator equations. 

The Controller Equation 

Substituting Eqs. (13), (14), (22) and (23) in Eq. 
(26), the state-estimation feedback controller equations 

dzJdt — [F + G W K W + GCg + DqBqI / s / z —DoB 0 y s (27) 

a = C 0 z (28) 


B 0 = -(PH/+ R wv )R v ~ J (22) 

D 0 = (I- y~ 2 PS)~ ] , p(PS) < y 2 (23) 

dP/dt = PF t + FP + G w R w G w r - (PH s r + R wv ) R V ~ J 
(PH s t + R wv ) t + P(y~ 2 Q X )P (24) 

where R v = E SW R W E S W T and must be positive definite and 
R wv = G W R W E S J. In Eq. (23) the spectrum p[PS] must 
me less than y 2 for D 0 to exist. The positive definite 
symmetric steady-state solution for P in Eq. (24) is 
obtained from the (2n x n) eigenvectors of the n stable 
eigenvalues of the estimator Hamiltonian matrix, 

(f - r kt r; 1 h 2 ) t (-//J/< // 2 + v- 2 // 7 ^, 1 //,) 
-G w KGl - RJi'fC ~(F - R r;'h 2 ) 

(25) 

which is dual to the state regulator Hamiltonian matrix 
inside [ ] in Eq.(12). If the eigenvectors are partitioned 
into two (n x n) matrices X and A, then P = X~ 2 A. The 
state estimate vector z is given by the Eq. (26). 

dzJdt = F z + G w w +Gu + D 0 B 0 (H2z -y2) (26) 

The complete duality relations between the state regulator 
problem and the state estimator problem are presented in 
Table 1. 

state regulator state estimator 
F F t 

G H s t 


are obtained. The standard LQG solution is obtained 
when y = 00 . The minimum value of y 2 for which a stable 
solution exists for P in Eqs. (22-24), provides the 
minimax control law that minimizes the maximum 
singular value of the matrix [y s T Qys 1/2 u T Q u 112 ], for the 
closed loop system. 

Unified Design Procedure 

In this design, the output y,i was chosen to be the 
linear combination of the trailing edge and leading edge 
accelerometer output (zte-0.4 zle), same as that used in 
the final classical design. One advantage of this choice 
was that the plant had no transmission zeros in the open 
right half plane. Usually in a frequency domain H-infinity 
design, the plant equations are augmented with weighting 
transfer functions. In this time domain formulation, the 
weights were chosen as constants. These weighting 
constants are chosen as inverse of the desired magnitude 
of the weighted quantities. The initial controller was 
designed with a large value of y 2 , using the plant Eqs. (1- 
3), at 225 psf, in air, at Mach 0.5, assuming G w = G. The 
block diagram for this unified design procedure is shown 
in Figure 11. The detailed design steps are described next. 

State -feedback Regulator Design 

Initially, the maximum output of the 
accelerometer sensors were of the order O.lg, (see figure 
3), and control surface maximum root-mean-square 
deflection was desired to be of the order 1 degrees. Thus, 
the initial values of the weighting matrices were chosen as 
follows: Q x = [ H/Q ys H s ] , Q ys = [100], and Q u = [1] . 
It was interesting to note that, instead of setting the cross 
weighting matrix Q xu = [0] as usual practice, the cross 
weighting matrix Q xu can be selected to place all state 



regulator poles beyond a certain distance a to the left of 
imaginary axis. This selection is accomplished by using 

Qxu = - aC, u (GjG u r l Qu (29) 

so that in Eq. (12), the eigenvalues of the diagonal matrix 
block F are off-set by 

GQu'Qmi T =~ al. (30) 

The control-weighting matrix Q u was 
subsequently reduced to 0.01 after a few design cycles to 
improve the regulator performance. This process of 
reducing Q u is equivalent to the state estimator loop- 
transfer recovery technique at the plant output. 



Fig. 11. Unified minimax control law design and 
evaluation procedure block diagram. 

State-estimator Design 

The state estimator was designed as a dual to the 
state-regulator with R w = 1, each diagonal elements of 
R v = 0.01, and R wv = [0], After a few design cycles the 
performance of the combined full-order controller was 
examined, and then R w was increased to 36. Since we 
also choose G w = G, this was equivalent to asymptotic 
state-regulator loop-transfer recovery at the plant input. 

4th Order Optimal Control law 

Subsequent solutions to the state-regulator and 
state-estimator were obtained with the same choice of 
weighting matrices but for decreasing value of y^, for 
which positive definite solutions for S and P could be 
obtained. Note that feasible solutions can be obtained for 


lower values of y^ up to y^> p (PS), below which the 
disturbance authority exceeded the control authority. The 
4th-order optimal control law was designed with y 2 =50 in 
order to obtain a low bandwidth controller. Fig. 12 shows 
the key singular value plots for analysis of multivariable 
stability margins to multiplicative and additive 
perturbation A at the plant input and output, with this 
minimax optimal control law, denoted as control law 4 . 
The minimum singular value o(I+KG ) is increased to 0.9 
at plant input and at plant output o(I+GK) is increased to 
0.5 from the corresponding values of 0.8 and 0.3 for 
control law 3, shown in Fig. 10. 


-(l+A)G-> -G(l+A)-> 



Fig. 12 Singular value plots for analysis of multivariable 
stability margins with minimax optimal control law 4, at 
225 psf, in air. 


Open Loop 225 psf Closed Loop, Law 2 225 psf 



Time, sec Time, sec 


Fig. 13 Open-loop and closed-loop responses due to step 
input 6te, with classical control laws 2, and 3, and 
minimax optimal control law 4, at 225 psf, in air. 

Fig. 13 shows the open-loop and closed-loop 
responses due to unit step input 6te of trailing edge 
control surface, at 225 psf, in air, at Mach 0.5, with initial 
control law 2, classical control law 3, and minimax 
optimal control law 4. The transient responses indicate 
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that the classical control law 3 provided better damping 
with lower control surface activity although the 
minimax control law 4 provided better robustness 
properties. This is the traditional trade-off between 
performance and robustness. The classical control law 3 
was implemented and tested in wind tunnel. These test 
results along with those using two optimal control laws 
designed by Waszak 13 are presented next. 

Flutter Suppression Test Results 

The performance and robustness of the final 
design was tested using the original full plant state-space 
equations and filters required for digital implementation. 
The 25 Hz antialiasing filters 157/(s+157) were added to 
the plant output. The washout filter 55/fs+5J and 
computational delay were added to the controller output 
equations. The 1/200 second computational delay was 
modeled by a (400-s)/(400+s) filter. Before the wind- 
tunnel test entry, the digital implementation was also 
numerically simulated. The numerical simulation block 
diagram of the control system using the final classical 
control law 3 is shown in Figure 14. This nonlinear 
simulation also included the effects of a dead-band 
present in the electro-hydraulic actuator. Application of 
the upper and lower spoiler for transonic flutter 
suppression with the same digital control law was also 
investigated using this simulation. 

The active flutter suppression control-law using classical 
design was successfully tested in air and in heavy gas 
medium at transonic speeds up to Mach 0.95. The tests in 
air indicated an increase in the flutter instability boundary 
from the open-loop dynamic pressure of 158 psf (Mach 
0.38) to the tunnel limit of 200 psf. A summary of flutter 
suppression test results in heavy gas is shown in Fig. 15. 
The solid line indicates the experimental flutter boundary, 
with the transonic dip at Mach 0.8. The tests at Mach 0.8 
indicated an increase in the flutter stability boundary from 
the open-loop dynamic pressure of 142 psf to the tunnel 
upper limit of 200 psf. A non-design plunge instability 
condition was also successfully suppressed. Classical 
control law 3 exhibited superior performance and was 
demonstrated to be stable with gain variation from 0.25 to 
7, and phase variation from -90 to 60 degrees. A non- 
design plunge instability condition was also successfully 
suppressed. Comparison of open-loop and closed-loop 
root mean square (RMS) responses of trailing edge 
accelerometer and control surfaces using the present 
classical control law 3 and two other control laws 
designed by Waszak 13 are shown in Figs. 16 and 17, 
respectively. These two control laws used upper and 
lower spoilers as control surface, for flutter suppression. 
Fig. 16 indicates that when the system is open loop stable, 
closing the loop actually reduces the response by 30%. 
Fig. 17 indicates that that the classical control law 3 
generally requires less control activity. All three control 


laws are comparable in performance, with control law 3 
exhibiting higher stability margins. 



Fig. 14 Numerical simulation block diagram of the control 
system digital implementation using the final classical 
control law 3 for flutter suppression. 


Closed Loop Test Points 



Fig. 15 Open-loop flutter boundary and closed-loop flutter 
suppression results from wind-tunnel tests in heavy gas. 



125psf 150psf 178psf 195psf 


Fig. 16 Open-loop and closed-loop RMS responses with 
classical control law 3 and two control laws employing 
spoilers, from wind-tunnel tests in heavy gas medium. 
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Conclusions 


Simple classical control laws, when properly 
designed based on physical principle, can successfully 
suppress transonic flutter and provide significant stability 
robustness in presence of shock and flow separation. 
Comparable robust optimal control laws can also be 
designed using a new generalized unified minimax 
formulation. Verification and improvement of the 
multivariable system stability robustness to unstructured 
perturbations at the plant, input and output were important 
steps in such a design process. Wind-tunnel tests in air 
and heavy gas indicated an increase in the transonic 
flutter dynamic pressure to the tunnel limit upper limit of 
200 psf. The control law robustness and performance 
predictions were verified in highly nonlinear flow 
conditions, gain and phase perturbations, and spoiler 


deployment. 



M= 0.63 M= 0.71 M= 0.77 M= 0.83 

125psf 150psf 178psf 195psf 


Fig. 17 RMS responses of the control surfaces with the 
classical control law 3 and two control laws employing 
spoilers, from wind-tunnel tests in heavy gas medium. 
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Appendix 

Reduced 4 th order state space equations in air at 225 psf. 
and Matlab script for unified minimax formulation and 
solution 


F = [ -1.6073 21.0010 0. 

-21.0010 -1.6073 0. 

0. 0. 0.7515 

0. 0. -25.1670 


0 . 

0 . 

25.1670 
0.7515 ] 


[G G w ] = [ -3.8259 
12.7130 
- 2.2202 
4.1351 


0.0597 
0.2720 
-0.1107 
-0.1745 ] 
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H d = [ -0.0517 -0.0132 -0.0668 0.0063 

-0.0542 -0.0090 -0.0482 -0.0016 

0. 0. 0. 0. 

0. 0. 0. 0. 

0. 0. 0. 0. 

10.3780 0.6369 7.9924 0.0671 

0.3897 -0.9373 -2.7625 0.9909] 

[E du E dw ] = [0.0440 0.0002 

0.0421 0.0004 

1.0000 0. 

50.0000 0. 

0. 0.0968 

0.5758 -0.0004 

0.0358 -0.0003 ] 

% Matlab script for unified formulation and solution 
% xdot = F x + Gw w + G u 
% yd = Hd x + Edw w + Edu u Design output 

% ys = H x + Esw w + Esu u Sensor output 

% 7 design output [zte zle dte ddte gust lift moment] 

% 

Qll = h'*qh*h+hd'*qhd*hd; 

Q22 = [q2]+[Edu'*qhd*Edu]; 

% pole placement using cross weight Q12 
alpha = 3.0; 

Q12=-alpha*g*((g'*g)\Q22)+[hd'*qhd*Edu]; 

WW = [Qll Q12 ; Q12' Q22]; 

% Generalized LQR Weights 
Ru = [gw*rw*gw'+nu*g*g']; 

Rv = [rv]+[Esw*rw*Esw']; 

Rwv = [gw*rw*Esw']; 

% generalized EST Weights 
VW = [Ru Rwv ; Rwv' Rv ]; 

% Generalized DESIGN with Edw=0 
[af,bf,cf,df]=lqg(f,g,h,Edu,WW,VW); 

% 

% STATE REGULATOR 

% Check if Q 1 1 is positive semi-definite and symmetric 
if any(eig(Qll) < -10*eps) I (norm(QlT-Qll.l)/norm 
(Q1 1,1) > eps) 

error('Qll must be symmetric positive semi-def), end 
% Check if Q22 is positive definite and symmetric 
if any(eig(Q22) <= 0) I (norm(Q22'-Q22,l)/norm(Q22,l) 
> eps) 

error('Q22 must be sym. positive definite'), end 

% 

% Construct Hamiltonian 

Hm=[(f-g*[Q22\Q12']) -g*[Q22\g']+gw*[rw\gw , ]/mu 
-Q1 1-Q12*[Q22\Q12'] -(f-g*[Q22\Q12'])']; 
[v,d]=eig(Hm); 

% 

% Sort eigen vector of neg eigenvalues 
d = diag(d); 

[d.index] = sort(real(d)); 
if (~( (d(n)<0) & (d(n+l)>0) )) 


error('Can"t order eigenvalues'), end 
% select vectors with negative eigenvalues 
chi = v(l:n,index(l:n)); 
lambda = v((n+l):(2*n),index(l:n)); 

S = real(lambda/chi); 

% 

% Positive feedback gain ku 

ku=-Q22\(g'*S+Q12'); 

kw=[rw\gw']*S/mu; 

% 

% closed loop plant f = fcl 
fcl=f+g*ku+gw*kw; 

X=lyap(fcl,eye(n)) % assume xo*xo'=I 

Xg=lyap(fcl,gw*gw'*36) 

W=0.5*trace(kw'*rw*kw*Xg) 

Jo=0.5*trace(S) 

Uu=trace(ku'*ku*Xg) 

Ycov=Hd*Xg*Xg'*Hd' 

% 

% H-inf ESTIMATOR DESIGN 
% define Hamiltonian Jam 
Ru = [gw*rw*gw'+nu*g*g']; 

J = [ (f-[Rwv/Rv]*h)' [-h'/Rv]*h+[hd7qhd]*hd/mu 
-Ru-[Rwv/Rv]*Rwv' -(f-[Rwv/Rv]*h) ]; 
[v,d]=eig(J); 

% 

% sort eigenvector of stable eigenvalues 
d = diag(d); [d, index] = sort(real(d)); 

% sort on real part of eigenvalues 
if (~( (d(n)<0) & (d(n+l)>0) )) 
error('Can"t order eigenvalues'), end 

% 

% select vectors with negative eigenvalues 

chi = v(l:n,index(l:n)); 

lambda = v((n+l):(2*n),index(l :n)); 

P = real(lambda/chi); 
ky = -(P*h'+Rwv)/Rv; 

% check positive definiteness of kd 

% Is this matrix (I - P*S/mu) nonsingular ? 

mumin=max(abs(eig(P*S))) 

if (mu > mumin), 

kd = inv(eye(n)-(P*S)/mu); , else 

error('spectrum(P*S) > mu , increase mu'), end 

% 

eve=eig(f+kd*ky*h); 

% controller structure 

Ao=f+g*ku+gw*kw+kd*ky*h; 

evc=eig(Ao); 

% H-infinity controller 
% Kop = [Ao -kd*ky -ku zeros(nc,ns)]; 

[fc, gc ,hc, ec] = feedback(f, g, h, e, Ao, -kd*ky ,-ku, 
zeros(nc, ns)); 

% 
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